//
// finFit.cc
//
// fit financial quotes distribution as
// a Gaussian or other bell shaped distribution whose mean is
// an exponential function of time (= constant interest rate)
// and whose width is a fixed fraction ot its mean
// BEWARE: fit convergence is extremely problematic
// 
// Requirements:
// - root v3.10
// - RooFit 1.04
// - finReadCsv.cc, finReadCsv.hh
//
// Usage:
// > root
// root> gSystem->SetIncludePath("-Wno-deprecated -Wno-overloaded-virtual");
// root> gSystem->Load("libRooFitCore.so");
// root> gSystem->Load("libRooFitModels.so");
// root> .L finReadCsv.cc+
// root> .x finFit.cc+
//

#include <stdlib.h>

#include "TCanvas.h"
#include "TH2F.h"

#include "RooFitCore/RooRealVar.hh"
#include "RooFitCore/RooGenericPdf.hh"
#include "RooFitCore/RooDataSet.hh"
#include "RooFitModels/RooGaussian.hh"

#include "finReadCsv.hh"

//
// fit options
//
// "m" = MIGRAD only, i.e. no MINOS
// "s" = estimate step size with HESSE before starting MIGRAD
// "h" = run HESSE after MIGRAD
// "e" = Perform extended MLL fit
// "0" = Run MIGRAD with strategy MINUIT 0 (faster, but no corr. matrix at end)
// 
// "q" = Switch off verbose mode
// "l" = Save log file with parameter values at each MINUIT step
// "v" = Show changed parameters at each MINUIT step
// "t" = Time fit
// "r" = Save fit output in RooFitResult object

// normalized Lorentzian	   
// L(x,x0,s) = 1/pi * 0.5*s / ( (x-x0)^2) + (0.5s)^2 )

// std::vector<quote_s> quotes; ---- not supported in root 3.10!

//--- maximum number of daily quotations
#define MAXQUOTES (6000)

//
// fit quotations exponential f. of time with a random spread
// proportional to the predicted quotation
//

int fitQuotes ()
{
  //
  // structure to store quotations data
  //
  quote_s qq[MAXQUOTES];
  vquotes_s vquotes(qq, MAXQUOTES);
  
  TString fname("../finData/VFINX.csv");
  // TString fname("../finData/VWESX.csv");
  
  if (readCsvQuotes(fname, vquotes)) {
    return EXIT_FAILURE;
  }

  //--- dataset variables
  RooRealVar quote("quote", "quote", 10, 160, "$");
  RooRealVar time("time", "time", 15, 40, "pyears");

  //--- parameters
  RooRealVar q0("q0", "q0", 3, 30, "$");
  q0.setVal(8);
  RooRealVar grate("grate", "grate", 1, 20, "%/pyears");
  grate.setVal(6);
  RooRealVar sigma("sigma", "sigma", 0.05, 0.8, "");
  sigma.setVal(0.2);

#if 1
  //
  // fit with generic PDF:
  // - 1 / ( (quote-qp)^2 + (gsigma/2)^2 )
  //   - with gsigma = qp*sigma
  //   - with qp = q0*exp(grate/100*time)
  // - Lorentzian errors
  // - Lorentzian sigma proportional to quotation
  // 

  // grate.setConstant();
  // sigma.setConstant();
  RooGenericPdf gp("gp", "quote vs. time",
		   "1/(pow(quote-q0*exp(grate/100*time),2)+"
		   "pow(sigma*q0*exp(grate/100*time)/2,2))",
		   RooArgSet(quote,q0,grate,time,sigma));
#endif

#if 0
  //
  // fit with generic PDF:
  // - 1/gsigma * exp(-1/2 * ((quote-qp)/gsigma)^2
  //   - with gsigma = qp*sigma
  //   - with qp = q0*exp(grate/100*time)
  // - Gaussian errors
  // - Gaussian sigma proportional to quotation
  // 
  RooGenericPdf gp("gp", "quote vs. time",
		   "1/sigma/q0/exp(grate/100*time)*exp(-0.5*"
		   "pow((quote/q0/exp(grate/100*time)-1)/sigma,2))",
		   RooArgSet(quote,q0,grate,time,sigma));
#endif

#if 0
  //
  // fit logarithm of quote rather than quote with generic PDF:
  // - 1/gsigma * exp(-1/2 * ((quote-qp)/gsigma)^2
  //   - with gsigma = qp*sigma
  //   - with qp = q0 + grate/100*time
  // - Gaussian errors
  // - Gaussian sigma proportional to log(quotation)
  // 
  RooGenericPdf gp("gp", "quote vs. time",
		   "1/sigma/(q0+grate/100*time)*exp(-0.5*"
		   "pow((quote/(q0+grate/100*time)-1)/sigma,2))",
		   RooArgSet(quote,q0,grate,time,sigma));
#endif

#if 0
  //
  // try to use RooFormulaVar:
  //   does not work when formula parameters are in dataset
  //
  RooFormulaVar mean("mean", "q0*exp(grate/100*time)",
		     RooArgList(q0,grate,time));
  RooFormulaVar gsigma("gsigma", "mean*sigma",
		       RooArgList(mean,sigma));
  RooGaussian gp("gp", "quote vs. time", quote, mean, gsigma);
#endif

  //--- dataset
  RooDataSet data("data", "data", RooArgSet(quote,time));

  //
  // fill dataset
  //
  for(unsigned int i=0; i < vquotes.size; i++) {
    quote.setVal(vquotes.q[i].quote);
    time.setVal(vquotes.q[i].fTime);
    data.add(RooArgSet(quote, time));
    // std::cout
    //   << vquotes.q[i].quote << " " << vquotes.q[i].fTime << std::endl;
  }
  
  gp.fitTo(data, "0q");

  TCanvas* canvas = new TCanvas("c1", "c1", 800, 800);
  canvas->Divide(1,2);

  TH2F* dh2 = time.createHistogram("quote vs time data", quote);
  data.fillHistogram(dh2,RooArgList(time,quote));
  canvas->cd(1);
  dh2->Draw("LEGO");

  TH2F* ph2 = time.createHistogram("quote vs time PDF", quote);
  gp.fillHistogram(ph2,RooArgList(time,quote));
  canvas->cd(2);
  ph2->Draw("SURF");

  return EXIT_SUCCESS;
}

//
// main
//

void finFit()
{
  fitQuotes();
}
